RBC-UKQCD Collaboration 



Proton lifetime bounds from chirally symmetric lattice QCD 

Y. Aoki, 1 P. Boyle, 2 P. Cooney, 2 L. Del Debbio, 2 R. Kenway, 2 CM. Maynard, 3 A. Soni, 4 and R. Tweedie 2 

(RBC-UKQCD Collaboration) 

1 RIKEN-BNL Research Center, Brookhaven National Laboratory, Upton, NY 11973, USA. 
2 SUPA, School of Physics, The University of Edinburgh, Edinburgh EH9 3JZ, UK 
3 EPCC, School of Physics, The University of Edinburgh, Edinburgh EH9 3JZ, UK 
' ''High Energy Theory Group, Brookhaven National Laboratory, Upton, NY 11973, USA 

(Dated: June 5, 2008) 

We present results for the matrix elements relevant for proton decay in Grand Unified Theories 
(GUTs). The calculation is performed at a fixed lattice spacing a -1 = f.73(3) GeV using 2+1 
flavors of domain wall fermions on lattices of size 16 3 x 32 and 24 3 x 64 with a fifth dimension 
\ of length 16. We use the indirect method which relies on an effective field theory description of 

. proton decay, where we need to estimate the low energy constants, a = —0.0112(25) GeV 3 and 

P = 0.0120(26) GeV 3 . We relate these low energy constants to the proton decay matrix elements 
£^ ' using leading order chiral perturbation theory. These can then be combined with experimental 

^ . bounds on the proton lifetime to bound parameters of individual GUTs. 
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& ; I. INTRODUCTION 



Proton decay is a distinctive experimental signature of Grand Unified Theories (GUTs). Decay experiments can 
test the predictions of these theories, and, even though direct nucleon decays have not been observed, the experimental 
lower bound on the decay rate has already ruled out the simplest minimal supersymmetric models [l[. One of the 
expected decay channels is N — > M + I, where N and M indicate respectively the nucleon and a pseudoscalar meson 
(K, 7r), while I is a lepton (e, fi, v e , v^). This decay is induced by supersymmetric particles or heavy gauge boson 
exchange, which can be integrated out to obtain an effective Lagrangian describing the low-energy dynamics in terms 
U^ ■ of the usual Standard Model fields. The lowest-dimensional operators that appear in this approach are (q c q){l c q) 
operators of dimension six; the transition amplitude for the decay is proportional to their hadronic matrix elements 
(M\(q c q)(l c q)\N) . A quantitative estimate of such hadronic matrix elements, which requires taming non-perturbativc 
QCD effects, is a key ingredient in probing the effects of higher-dimensional operators in GUTs models at current 
i experiments. 

OO ' Several determinations of the hadronic matrix elements have been performed in the past @, [E 0, [E @, 0, H, IE EE 
El El, EE USE!, usm S either QCD bound state phenomenological models, or lattice QCD. Due to recent progress in 
simulating dynamical fermions, lattice QCD has become a quantitative method to compute hadronic matrix elements 
from first principles with controlled systematic uncertainties. The matrix elements relevant for nucleon decay can be 
■ extracted from three-point correlators computed on the lattice. This is the so-called direct method, which requires an 
expensive computer simulation. However, the same matrix elements can also be computed, at the cost of introducing 
difficult to estimate systematica, using the chiral lagrangian describing proton decay [16|: in this case they are 
expressed as functions of the low-energy constants (LECs) that appear in the chiral lagrangian. These LECs can be 
computed from two-point lattice correlators at a lesser computational cost. However, such an indirect determination 
of the matrix elements depends on the accuracy of chiral perturbation theory, and therefore is affected by an additional 
source of systematic error. 

In this work, we present a new determination of the matrix elements that are relevant for nucleon decay based on 
the indirect method, using dynamical Domain Wall Fcrmion (DWF) configurations with 2 + 1 flavours. Our results 
extend the ones obtained for 2 flavours of dynamical DWF in Ref. [l^|. Systematic errors are greatly reduced by 
simulating at light quark masses, and using non-perturbative renormalization. Note that the exponentially suppressed 
chiral symmetry breaking of DWF greatly simplifies the mixing of operators under renormalization, which improves 
the precision of the final result. The correct number of flavours gives confidence in setting the scale, a large source of 
uncertainty in some early lattice determinations. 

The chiral perturbation theory results used for this work are summarized in Section [Hi which also sets the notation 
used throughout the paper. The chiral Lagrangian for nucleon decay involves two LECs, which are obtained by 
extrapolating to the chiral limit the outcome of numerical simulations performed at light quark masses. A direct 
measurement of the hadronic matrix element using lattice three-point functions, which relies much less on the validity 
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of chiral perturbation theory, is in progress, and is deferred to a subsequent publication. 

Details of our lattice simulations are reported in Section IIIII Our results are obtained from gauge configurations 
with volumes of 16 3 x 32 and 24 3 x 64. Both have a fifth dimension of size L s = 16 and use the Iwasaki gauge action. 
The lattice spacing is a w 0.114 fm corresponding to physical volumes of (1.8 fm) 3 , and (2.7 fm) 3 respectively. Being 
able to compare two different physical volumes enables us to estimate the finite volume effects. Results for meson 
spectroscopy and topology have already been presented in Ref. for the smaller lattice and in Ref. [HI for the 
larger lattice. We refer to these publication for details of calculations involving the lattices which are used in this 
paper. The range of fermion masses simulated for this work yields a ratio of the pscudoscalar to vector meson masses 
in the range 0.378 < mps/my < 0.615. Working at fixed lattice spacing, we are not able to present a continuum 
extrapolation for our final result. Nonetheless, it should be noted that DWFs are automatically 0(a) improved [4l| . 
and is therefore expected to have a good scaling behavior. 

Section HVl presents our results for the non-perturbative renormalization (NPR) of the three-quark operators, using 
the RI-MOM scheme. We compute the renormalization mixing matrix and perform the matching required to obtain 
the renormalized operators in the MS-scheme. 

In the last Section of the paper, we combine the lattice amplitudes with the renormalization factors to compute 
the phcnomcnologically relevant matrix elements in the MS scheme. We discuss the error budget in detail including 
estimates of the systematic error due to the chiral extrapolation, the renormalization, the finite volume and the choice 
of method to set the lattice spacing, and the foreseeable improvements on the current estimate. 



II. CHIRAL LAGRANGIAN FOR PROTON DECAY 



Integrating out the heavy GUT particles yields the low energy effective Lagrangian describing nuclcon decay written 
in terms of the QCD fundamental fields: 



(i) 



d=1.2 i=l 



rf=l,2 t=l 



where d denotes the generation of the lepton produced in the decay, and i is a label for the dimension-six operators 
containing three quark and one lepton field that describe nucleon decay. C^' and are Wilson coefficients. The 
full list of operators was identified on symmetry grounds in Refs. 0, H0,[ll|; their matrix elements between 

hadronic states determine the decay amplitude. For instance, the matrix elements that are relevant for the process 
where a proton decays into a pion are: 



(7r(p)\e abc (u al CP R ,Ld b )P L u c \p(k, S )) = P L 



R/LL, 



{q*)-iiW*l LL {q 2 ) u{k,s), 



(2) 



where a, b, c are colour indices, C is the charge-conjugation operator, and Pr.l = 1= 2 75 are the right- and left-handed 
projectors. The non perturbativc dynamical effects are captured by the two form factors that appear on the RHS 
of Eq. @, while q (the momentum carried by the electron) is the momentum transfer. It is convenient to introduce 
here a generic notation for three-quark operators with an arbitrary spin structure: 



O rr '(x,t) = e abc [u a (x,t)(CT)d b {x,t)] T'u c (x,t). 



(3) 



where T and V are elements of the Clifford algebra in four-dimensional Euclidean spacetime, and we have omitted 
spinor indices. We use the notation S = 1, P = 75, V = 7 M , A = 7^75, T = = ^{7^,7^}, R = Pr, and 
L = Pl. Further operators with this structure appear when computing the nucleon mass, and upon renormalization, 
as discussed in Sects. IIIII and IIVI 

Following the notation in Refs. [l3l . [l6j . the chiral Lagrangian describing baryon-meson dynamics is written in 
terms of a pseudoscalar meson (octet) field: 
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and a spinor baryon (octet) field: 
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(5) 



At lowest order in powers of momentum, and in Euclidean space-time, the chirally symmetric Lagrangian is written 
as: 

Co =^-Tr(d^)(d^) + TcBfadp + M B )B 

o 



(£> + F)TrB lfll5 - ^d^] B, 



where the unitary matrices £ and £ are defined as 

£ = cxp 

Introducing a diagonal quark mass matrix: 



2i4 
T 



'i<f> 

£ = exp — 



(6) 



(T) 



M 



"id 



(8) 



the symmetry-breaking part of the chiral Lagrangian becomes: 

Ci = - u 3 Tr {pM + MY) - a{YxB {$M$ + £M£) S - a 2 TrBB ((*Mg + £M£) 
- 6iTrS 75 (^M^ - £M£) S - b 2 TrB~f 5 B - £M£) . 



(9) 



The low-energy constants that appear in the chiral Lagrangian are extracted from phenomcnological analyses. In 
particular, following the notation in Ref. [Hj], / is the pion decay constant in the chiral limit, 130(5) MeV [22|. The 
combination F + D yields the nucleon axial charge, qa = 1.2695(29) [22| . while the combination F — D is related 
to the ratio of the zero-momentum form factors for semilcptonic hyperon decay, <7i//i [23]. Together these give, 
F = 0.47(1) and D = 0.80(1). a\ and a 2 are symmetry-breaking parameters, but their values are not required in this 
work. The parameters 61,62 are not precisely determined and are an extra source of systematic error. 

The transformation properties under SU(3) l x SU(3),r of the three-quark operators in Eq.Q]determine the expression 
of the baryon-number violating operators in the chiral Lagrangian. The latter appear in the Lagrangian with two 
new low-energy constants a and /3 [l6| : 



d=l 



+ C d 2) e dR TvT^B R ^ + 

+cf ) e dR TrFgB R g + cf \ iL ^T"iB L ^ 
2 

{ C d ] h L Tr^B L ^ - v dL TrT'iB L e] 



e dL TrFZB L £ - v dL TiP^BU 



d=l 



C d i} e dR Tv^B R i + C\ 



e dL Ti^B L ^ - VdL TxPiB L S} 



+C^e dR Tr^B R ^ + cf ] v dL Trf" £B L £}} + h.c. 



(10) 
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The matrices J 7 , J 7 ', J 7 , J 7 ', and J-" are projectors in flavour space; their explicit expressions are: 



J 7 = 








Eqs . 1 1 01 and 1111 show that the low-energy constants a and (3 determine the matrix elements: 

(0\O RL \p(k,s)) =a P L u(k,s) 
(0\O LL \p(k,s)) =/3P L u(k,s) 



(11) 



(12) 
(13) 



where u(k, s) is the spinor associated with a proton of momentum k and spin projection s. The phase definition is 
fixed such that a and (3 are real and a < 0. As we will later describe, we observe a + (3 ~ 0, which is expected because 
of the relation, 



{a + 13) u(k,s) 



e abc (u la Cd b ) l5 u c \p{k,s)), 



(14) 



which vanishes in the non-relativistic limit and is known to be quite small even at small quark masses [24j . 

Using chiral perturbation theory to compute the matrix element in Eq. ([2|) yields for the N — ► n transition [l3L [l6[ : 



(ir \O RL \p(k,s)) = aP L u(k,s) 

— aPi,iq l u(k., s) 
(ir \O LL \p(k, S )) = (3P L u(k,s 

— (3PLi4u{k,s) 
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(15) 



(16) 



where q is the four-momentum of the outgoing lepton. In the limit where q 2 <C m 2 N and b\m u <C m^, these expressions 
simplify to: 



(n \O RL \p(k,s)) 

(Ao LL \ P {k,s)) 



uPLuik, s) 
0P L u(k, s) 



1 D + F 
1 D + F' 



V2f V2f 



+ 0(m 2 /m 2 N ) 7 
+ 0(mf/m 2 N ), 



(17) 
(18) 



-q" = mf is the on-shell condition for the outgoing lepton. The equations above relate the proton decay 
matrix elements to the low-energy constants a and /?; note that, in order to reconstruct the matrix elements on the 
Ihs of Eqs. im and [T8l using the indirect method, the combination F + D and the pion decay constant, /, are also 
required. 



III. LATTICE SIMULATIONS 
A. Dataset description 

The analysis was performed on 2+1 flavor DWF ensembles with volumes of 16 3 x 32 and 24 3 x 64 generated using 
the Iwasaki gauge action with (3 = 2.13 and the domain wall fcrmion quark action with L s = 16. At each volume we 
generated sets of configurations with a light isodoublet with masses am u d = 0.005 (24 3 x 64 only), 0.01, 0.02 or 0.03 
and a fixed approximate strange quark mass, am s = 0.04. The ensembles, described in [25| and [l8j . have a fixed 
inverse lattice spacing of a -1 = 1.73(3) GeV and were generated with the RHMC algorithm with a trajectory length 
of t = 1. These same datasets were used to calculate gA, [2(|. The configurations used for both the non-perturbative 
renormalisation and the matrix element calculation are shown in Table |TJ 

For each of the seven ensembles matrix elements were calculated using correlation functions composed of valence 
quarks with masses equal to the light quark mass in the sea. To improve statistics, correlators were oversampled and 
averaged into bins whose size depended on the Monte Carlo time separation between measurements. The binning was 
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V x L s 



am ud 



Matrix Elements 

iVtraj A iV c f g A sr c A^bir 



NPR 

iVtraj A iV c f g iV src iVbir, 



0.01 

16 3 x 32 x 16 0.02 
0.03 



500-4000 10 175 4 8 
500-4000 10 175 4 8 
500-7580 10 177 2 8 



1000-4000 40 75 4 1 
1000-4000 40 75 4 1 
1000-4000 40 75 4 1 



0.005 

24 3 x 64 x 16 0.01 
0.02 
0.03 



900-4500 10 90 2 8 

10 90 2 x 2 8 

1500-3860 10 59 2 8 

40 59 2 2 

1800-3600 10 45 2 8 

40 45 2 2 

1020-3060 20 51 1 2 

40 51 1 1 



TABLE I: RHMC 2+1 flavour datasets used for the non-perturbative renormalization and matrix element calculation. V is 
the space-time volume of the lattice, L B is the extent of the fifth dimension, am u d is the up sea quark mass (the strange sea 
quark mass is kept fixed at 0.04), iV tr aj is the lowest to highest trajectories analysed with matrix elements calculated every A 
trajectories, JV c f g is the number of configurations, iV S r C is the number of quark propagators solved with different source locations 
and iVbin is the bin size. 24 3 x 64 x 16 data were generated for the non-perturbative renormalization calculation, however, it 
was only used as a check for finite volume errors and so does not appear here. For the 24 3 x 64 x 16 matrix element data, 
there are two independent runs for each of the sea quark masses. These independent runs used different smearings, A, source 
locations and iV S rc- 

consistent with the integrated auto-correlation length for the pscudoscalar meson correlators at the time separation 
typically used, which was calculated to be of order 50 trajectories. Multiple sources per configuration and several 
different types of smearing have also been used to improve the signal. As well as local sources (L), we employ gauge- 
invariant Gaussian smearing with two different smearing radii (G and G*) and gauge fixed hydrogen-like wavefunction 
smearing (if ) . One or both of the propagators used to construct the two-point correlators for mesons may be smeared 
while for baryons, one, two or all three propagators may be smeared. We adopt the same convention used in Ref. [27j | 
for naming the smeared two-point functions. 

The chiral limit is defined as the value of arrif such that arrif + am res = 0, where am res = 0.00315(2) is the residual 
quark mass, estimated in Refs. (2f| H3, H3- The lattice scale is determined from a combination of the Vt~ baryon 
mass and the pscudoscalar kaon and pion masses, yielding a value a -1 = 1.729(28) GeV (see Ref. [l8[ for details). 

The parameters used in the simulations correspond to a pseudoscalar meson mass ranging from 331 MeV to 671 
MeV. The renormalization constant for the axial current, Za = 0.7162(2), which we will use in the non-perturbative 
renormalization of the nucleon decay operators, was obtained from a hadronic matrix element of the conserved DWF 
axial current in Ref. [28|. 

B. Nucleon mass and amplitude 

Starting from the correlator of two operators, TlV2 and Ts,Ti 1 we can define the scalar two-point function: 



/r 1 r 2 ,r 3 r 4 (i) = J2 Tr 



iivrw- ^7n r 3T4^ N ( 1 + 74 



(O ll ' l2 (x, t)0 " 4 (0)) 



(19) 



Using the notation introduced so far, O f (x,t) is the usual local proton interpolating operator: 

O ps (x, t) = e abc [u aT (x, t)C l5 d b (x, t)] u c (x, t), (20) 
and the large time exponential fall off of the correlator fps.ps is dictated by the nucleon mass: 

fps,Ps(t) = ^~ amNt G% + ---: (21) 
where Gjy is the overlap of the proton interpolating field to the normalized proton state: 

(0|O PS (0, 0)|p(k, s)) = G N u(k, s). (22) 




FIG. 1: (a) is an effective mass plot and (b) is an effective amplitude (Eq |26|) plot for the nucleon. Both are calculated on the 
24 3 x 64 dataset with am u = 0.01. The different colours in the effective mass plot correspond to different smearings. Datasets 
are labelled with the smearing (i.e. LL). Those datasets labelled with a 2 use the operator /a 4 s,a 4 s(£), the rest use fps,ps(t)- 
(c) is a linear extrapolation of the ground state mass to the chiral limit. 



The nucleon mass is obtained from the two-point functions fps,ps(t) an d fA i s,A 4 ,s(t)i each of them being computed for 
several smearing combinations. Firstly, for each two-point function and for each smearing combination, we calculate 
the effective mass: 



m N . cS (t) = log 



/(*) 



(23) 



where f(t) indicates the two-point function in any one of the channels used for the analysis. Results for the effective 
mass computed from both two-point functions, and for two different smearing combinations, are reported in Fig. [1] 
The agreement between the different channels within the error bars is clearly seen in the first plot on the left for the 
24 3 x 64 data with am u = 0.01. The effective mass can be fitted to the same constant m for each channel; correlations 
between different time-slices are taken into account by minimizing a correlated x 2 '- 



<n)2 



(m (rl) ) 



E 

t.t' 



» 



C 



(n)-l 



l N.c 



» 



(24) 



where G tt , is the covariance matrix: 



C 



(n) _ 



boot 



N hoot 

E 



(n,m) 
l AT,cff 



(*) 



An) 
l N.c 



«(*)> 



l N,cfi 



>) 
l N,cS 



(*')> 



(25) 



the index n represents a bootstrap resampling of the original data, and the index m represents a bootstrap resampling 
of the n th bootsamplc. A^oot is the number of bootstrap samples, rn^'^g (t) is the effective mass determined from the 



..th 



resampling of the n th bootsample and (m^ cS (t)) is the average of the effective mass over the m th resampling of 
the n th bootsample. 

All channels display a plateau for the effective mass, and the limiting values are compatible within the statistical 
errors. The smeared propagators reach the limiting value earlier, as expected, thus yielding a longer plateau for the fit 
to be performed. In order to increase the precision of the fit, all channels were fitted simultaneously to a single constant 
to; correlations between different channels are also taken into account in the construction of the covariance matrix. 
Note that with this fitting procedure several channels are fitted simultaneously without adding extra parameters, and 
the minimization in the one-dimensional parameter space can be performed analytically. 

The fit to the amplitude, Gat, is subsequently performed by defining an effective amplitude: 



G 2 



1 



N, 



s(t) = -fps,ps(t)exp{mt), 



(26) 



where m is the nucleon mass obtained in the fit described above, and we used the same notation as in Eq. [23] to 
denote the two-point function. Gat is then obtained from a fit to a constant by minimizing a fully correlated x 2 , as 



(a) 




0.8 ■ 



I 




(b) 



FIG. 2: Effective mass plots for the nucleon on the am u — 0.005, V = 24 3 x 64 ensemble. The different colours in the effective 
mass plot correspond to different smearings. Datasets are labelled with the smearing (i.e. LL). Those datasets labelled with 
a 2 use the operator fA i s,A i s(t), the rest use fps,ps(t). (a) shows the fit before scaling the errors and (b) shows the fit after 
rescaling the errors 



discussed above for the nucleon mass. Results for Gjv.cffW are displayed in Fig.QJb), where a long plateau is clearly 
visible. 

For all the fits presented here, the results of the minimization procedure are stable with respect to sensible variations 
of the fit range. All the correlators, smearings, and fit ranges are summarized in Table [Til Variations of the fitted 
parameters remain within their statistical error as the bounds of the fitting range arc shifted by ±1 timcslicc. 

For the case of the nucleon mass on the am u = 0.005, V = 24 3 x 64 ensemble (the ensemble with the lightest 
valence quark mass), there was some difficulty in judging exactly where the plateau for the effective mass started. 
Fitting to different time ranges gave incompatible results. To account for this, we performed a fit over a large time 
range, spanning the multiple potential plateaux. The incompatibility of the data was reflected in a poor value of x 2 
per degree of freedom (d.o.f ) of 4.3. In order to deal with this we rescaled the errors on all the points in the effective 
mass plot by y / x 2 /d.o.f and performed a second fit to this rescaled data. This gave a x 2 /d.o.f of 1, as expected, and 
a fitted mass compatible with the best fit value from before, but with a larger error. The fits to the effective mass on 
this ensemble before and after rescaling are shown in Fig. [2j 

The nucleon mass and amplitude are extrapolated linearly to the chiral limit. The result of the extrapolation for 
the nucleon mass on the 24 3 x 64 dataset is displayed in Fig. QJ c) . The results for the nucleon masses obtained from 
the fits are summarized in Tabic IIIII 



C. Low— energy constants 



As discussed in the previous section, the low-energy parameters a and (3 appearing in the chiral Lagrangian can 
be calculated at leading order through the proton to vacuum matrix elements: 

(0\O RL \p(k, a)) = aP L u(k, s), (0\O LL \p(k, s)) = (3P L u(k, s), (27) 



-{0\O LH \p(k,s)) =aP R u(k lS ), -(0\O HH \p(k, S )) = /3P R u(k,s), (28) 

where Eq. [28] is obtained from Eq. [27] by parity transformation. The low-energy constants are obtained from the 
asymptotic behaviour of ratios of two-point functions for large Euclidean time t: 

R a (t) = 2G N { RL ' Ps{ f). -> a, (29) 

R (t) = 2G N { LL ' Psi fl ->/?. (30) 
Jps.psyi) 
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24 3 x 64 
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Smearing 


o rr 


mjv 


Gjv 


a 





Smearing 


o rr 


mjv 


Gjv 


Q 





0.005 














LL 


QPS 


- 


9-12 


5-8 


5-9 
















HL 


QPS 


6-12 


_ 


4-10 


4-9 
















HL 


QAiS 


6-12 


_ 


_ 


- 
















G*L 


QPS 


6-12 


_ 


_ 


- 
















G*G* 


QPS 


6-12 


- 


- 


- 


0.01 


LL 


o PS 


9-12 


9-12 


5-8 


5-8 


LL 


QPS 


9-12 


9-12 


7-11 


5-10 




LL 


QAiS 


9-12 


_ 




- 


LL 


QAiS 


9-12 


_ 


_ 


- 




GL 


qPS 


9-12 


_ 


3-8 


3-8 


GL 


QPS 


8-12 


- 


7-11 


4-10 




GL 


QAiS 


9-12 


_ 






GL 


QAiS 


8-12 


- 


_ 


- 
















G*L 


qPS 


7-12 


- 


_ 


- 


0.02 


LL 


qPS 


9-12 


9-12 


6-12 


5-12 


LL 


qPS 


9-10 


9-12 


7-11 


5-10 




LL 


QAiS 


9-12 


_ 






LL 


QAiS 


9-10 


_ 


- 


- 




GL 


QPS 


8-12 


_ 


6-10 4-14 


HL 


QPS 


9-11 


_ 


7-11 


7-10 




GL 


QAiS 


8-12 


_ 






HL 


QAiS 


9-11 


_ 


_ 


- 




GG 


QPS 


8-11 


_ 






G*L 


QPS 


9-11 


_ 


_ 


- 




GG 


QAiS 


8-11 


- 


















0.03 


LL 


QPS 


10-12 


9-12 


7-11 


6-11 


LL 


QPS 


10-12 


9-12 


9-13 


9-11 




LL 


QAiS 


10-12 








LL 


QAiS 


10-12 






- 




GL 


QPS 


8-12 




7-11 


8-12 


HL 


QPS 


9-12 




9-13 


9-11 




GL 


QAiS 


8-12 








HL 


qA 4 S 


9-12 










GG 


QPS 


8-12 








G*L 


QPS 


8-12 










GG 


QAiS 


8-12 





















TABLE II: Smearings, operators and fit ranges used for the calculation of nucleon masses, nucleon amplitudes and matrix 
elements. 



A typical plateau obtained for R a is shown in Fig. [3] Two different smearing combinations were used in the analysis, 
which correspond respectively to a local and a smeared interpolating field (x, t) for the nucleon. They both yield 
consistent results, as shown in the plot. The values of the low-energy constants were obtained by fitting the data to 
a constant, for each value of the quark masses. As for the spectrum, x 2 is always defined taking into account the 
correlation between different time-slices. The results obtained from the fits are given in Table [LTT1 

The data points are extrapolated linearly to the chiral limit, as shown in Fig. [4] The data in the mass range studied 
in this work appear to be consistent with a linear behaviour, leading to a good fit for the chiral extrapolation. An 
uncorrelated \ 2 is used in this case, since the points at different values of the quark mass are produced by independent 
runs. 



IV. NON— PERTURB ATIVE RENORMALIZATION 
A. RI-MOM mixing matrix 

For the non-pcrturbative renormalisation of the proton decay matrix elements we employ the non-perturbative, 
MOM-scheme, renormalisation technique of the Rome-Southampton group [jpl as used by [Hj|,[l9[ and [3(j. The 
application of this technique to proton decay matrix elements is outlined in [151 ] which we briefly summarise. 

The operators, O rr , can be classified according to their symmetry properties under parity (V) and the so-called 
switching transformation (S). The result of such classification is summarized in Table Hvl In the presence of chiral 
symmetry breaking, operators that belong to the same sector mix under renormalisation. Concentrating on the iS _ 
sectors, the renormaliscd operators in the parity basis arc defined as: 

Of en = z£EO? aH , A,B = {SS,PP,AA}, (31) 




(b) 
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B 

Pi 0.005 



• LL 

• HL 



0.002 
16 



D.I II IX 



a 

Qi 0.005 



• LL 

• HL 



IIS 



t 

(d) 



FIG. 3: The ratio R a in Eq.[30]for the 24 3 x 64 dataset with am u = 0.005, 0.01, 0.02 and 0.03 respectively. The different colours 
correspond to different source smearing. Horizontal lines show the fit to the plateau. 
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FIG. 4: Linear chiral extrapolation for the ratios R a (a) and Rp (b) for the 24 3 x 64 dataset. 
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V x L s am u d/am s 


amjv a?ot a 3 (3 


0.03/0.04 
16 3 x 32 x 16 0.02/0.04 
0.01/0.04 
chiral 


0.908(6) -0.00695(19) 0.00719(21) 
0.819(8) -0.00605(31) 0.00606(30) 
0.722(19) -0.00478(43) 0.00511(47) 
-0.00349(64) 0.00369(63) 


0.03/0.04 
24 3 x 64 x 16 0.02/0.04 
0.01/0.04 
0.005/0.04 
chiral 


0.892(10) -0.00689(33) 0.00621(38) 
0.805(12) -0.00571(32) 0.00598(38) 
0.720(10) -0.00508(29) 0.00486(28) 
0.671(5) -0.00397(18) 0.00400(22) 
-0.00326(27) 0.00348(32) 



TABLE III: Results from fits described in this paper. The nucleon masses the LECs a and j3 are reported as a function of the 
quark masses, for both lattices used in this study. The results of linear chiral extrapolations are also reported in the last line 
of each column. All the results are given in units of the lattice spacing a ~ 0.12 fm. 





s~ 


5+ 




SS, PP, AA 


VV, TT 


■P+ 


SP, PS, -AV 


-VA, TT 



TABLE IV: Classification of the O rr operators according to their transformation properties under parity and switching. 



where A and B label the possible choices for IT' and Z$B is a 3 x 3 mixing matrix. The same mixing matrix 
renormalizes the operators in the sector V~ and V . The chirality basis, which contains the operators of interest for 
the nucleon decay matrix elements, consists of 

LL = ^(SS + PP)-^(SP + PS) (32) 

RL = -(SS - PP)--(SP- PS) (33) 

A(LV) = l -AA- l -(-AV), (34) 

hence, the mixing matrix in the chirality basis, 2nd, and in the parity basis are related via: 

Znd = 1~Znd1~~ 1 , (35) 

where 

' 1/4 1/4 

T= | 1/4 -1/4 ] . (36) 
1/2 

The mixing matrix in the parity basis is computed from the non-perturbative amputated three-quark vertex function 
of the operators in the S~ sector as a function of external leg momentum p after gauge fixing to Landau gauge. The 
number of configurations used in the non-perturbative renormalisation is given in Table [J The vertex function is 
defined as the amputated Fourier transform of the correlator of O a with three quark spinors: 

G a (p 2 Uc a0 - (S = e Q ' 6 ' c '(cr) Q , /3 T^,(g^ Q (p)Q^(p)Qy7(p))' ( 37 ) 

where 

Q&(P) = (SaS&r'S&lip), Si, a a {p) = [ dxe-*-* Si,l{x), (38) 



S{x) is the quark propagator, and T,V are the matrices that appear in O a . 
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FIG. 5: The mixing matrix M in Eq. (|39[1 in the chirality basis, IT' = {LL, RL, A(LV)}, as a function of the lattice momentum 
squared for the 16 3 x 32 lattices with am u = 0.01,0.02 and 0.03 (from left to right, respectively). The off-diagonal mixing 
between operators is highly suppressed. It is worthwhile to note that the mass dependence of the mixing matrix is very mild. 
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Introducing the mixing matrix: 

M AB (p 2 ) = Q A hc al3 7<5 ■ P B bc j3a &1 , (39) 

the renormalisation condition in the RI-MOM scheme reads: 

Z q z l 2 Z B %M CA = S BA (40) 

where Z q is the quark wavefunction renormalisation; a, b, c are colour indices and a, (3, 7 and S are spin indices 
associated with T, I~" respectively. The projection operators, 

' C (75C- 1 )^75 7 (42) 
> bc (j 5 ^C-y a ( l5 ^) s \ (43) 

are chosen such that the renormalisation condition in Eq. HOI is satisfied in the free field case: Z q = 1, Z^£ = 8 . 

Fig. [5] shows the mixing matrix, M AB , in the chirality basis as a function of external leg momentum. The set of 
momenta used to calculate the mixing matrix is defined by 

/2tt 2tt 2tt 2tt \ 
P = -j-n x , —n y , —n z , —n t (44) 

Ijy J-iz '-'t J 

where L x = L y = L z is the spatial size of the lattice and L t is the time extent. Combinations of (n x , n y , n Zl n t ) such 
that —2 < n x , n y , n z < 2 and ~4 < n t < 4 are chosen and then averaged into equal p 2 values. 

Operator mixing is induced by chiral symmetry breaking. The extent to which chiral symmetry is broken in the 
domain wall action is parameterised by the residual mass, am rcs , and the induced mixing is expected to be suppressed 
by a factor (am res ) 2 [3l|. It may be seen from Fig. [5] that, in the window of momenta for which contributions from 
both hadronic effects (low momenta) and contributions from discretisation effects (high momenta) are small, the 
chiral symmetry afforded by the domain wall fermions suppresses the mixing between different chirality operators and 
results in a mixing matrix which is essentially diagonal. Th is g reatly simplifies the calculation of the proton decay 
matrix elements compared to, for example, Wilson fermions [1 31 ] - 

The matrix Znd can be obtained from the relation M = Zq^ 2 Z^, as shown in Eq. 0D] which requires Z q to be 
computed. Instead, we remove the Z q dependence, and exploit the accurate determination of Za = 0.7162(2) at the 
chiral limit, which was computed from ratios of hadronic matrix elements in Ref. [T^], together with the average of 
the amputated local axial vector and vector bilinear currents, which allows the evaluation of the factor A A = Z q /ZA- 
Fig. [B]shows the average and difference of the amputated local axial vector and vector bilinear currents. The non-zero 
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difference may be taken as a measure of the systematic error of the renormalisation constant arising from the closing 
of the window where the RT-MOM NPR can be safely applied. It may be observed that for (ap) 2 > 1.7 there is < 1% 
effect, which is enhanced to 2% by extrapolation of {ap) 2 — > 0. 

The product (A A ) 3 ^ 2 M~ x yields Z^d/Z Z / 2 for each value of the sea quark mass, without having to deal directly 
with Z q . At finite lattice spacing, Za§ only scale dependence is due to the discretisation error, which starts at 
0(a 2 p 2 ). Finally the rotation to the chirality basis and a linear chiral extrapolation are performed, the latter may be 
done very precisely, as the mass dependence is extremely mild, as shown in Fig. \S\ As an example, the p 2 dependence 
of the LL element of the matrix Z^n/Z^ 2 is displayed in Fig. [6l 



B. Scheme matching and RG running 

In order to relate the lattice, MOM-scheme, matrix elements at scale p to those in the MS, NDR scheme at some 
scale y, we compute the factor 



Z MOM (p)' 



ND 



(P), 



(45) 



where Z MS (p) / Z MOM (p) is the matching factor from MS scheme to MOM scheme at a scale p calculated using 
continuum perturbation theory, and U MS (fj,;p) is the renormalization group evolution factor from scale p to /i in the 
MS scheme. The matching factor has been computed in Ref. [l5j : 



Z MS _ a. 

Z MOM - 1 + 4?r 



433 
180 



1123, „ _/587 317, „ 

^r ln2+ Hi8o-w ln2 



(46) 



where £ = as we work in Landau gauge. The MS evolution factor reads 



U^{^p) 



Otsijp) 



7o/2/3o 



2 38 
A) = II-3W/, ft = 102 - —Nj, 



7o 



-4, 71 



14 



9 



iV> - 4A 



in 



(47) 
(48) 
(49) 



where the anomalous dimension of the nuclcon decay operator has been calculated up to two loops in MS, NDR 
scheme [32| and A = 0,-10/3 for LL,RL operators respectively. The value of a s (p) is obtained by integrating 
numerically the four-loop [3 function of Ref. [33|], starting from a s (Mz) = 0.1176(2) [231, an( i matching the value of 
a s across the b, and c thresholds. 

The MS renormalisation factor, Eq. 3S1 at a fixed scale /j, = 1/a is plotted in Fig. [S]as a function of the square of 
the scale at which the lattice, MOM-scheme, renormalisation calculation was performed. The remaining momentum 
dependence, due to 0(a 2 p 2 ) discretisation errors, is removed by performing a linear extrapolation in (ap) 2 to (ap) 2 = 0, 
which is also shown in Fig. [SJ This extraploation is performed over the range 1.7 < (ap) 2 < 2.5 where the non- 
perturbative effect, estimated at 2%, is expected to be small. 

Together with the value of Za from the hadronic matrix element ratio and using Eq. [47] to run from fj, = 1/a to 
[i = 2 GeV we obtain: 



U 
U 



MS^latt 



MS^latt 



(H = 2 GcV)ll 
(H = 2 GeV) RL 



0.662 ±0.010 
0.665 ±0.008 



where the error is statistical. 



V. DISCUSSION 



The errors on all quantities so far have been purely statistical. From the results in Table [TTT1 we can see that for 
this matrix element and for the statistics available, there are no significant finite volume effects as the results on both 
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FIG. 6: The figure on the left shows the average and the difference of the amputated local axial vector and vector bilinear 
currents as a function of (ap) 2 . In the figure on the right, the black points show the MOM-scheme renormalisation factor in the 
chiral limit for the O ll operator normalised by the axial current renormalisation factor Za as a function of the renormalisation 
scale (ap) 2 . The red points show the renormalization factor in the MS scheme at a scale (i = 1/a asa function of the matching 
scale. The red line shows the linear extrapolation in (ap) 2 , where the blue points are those included in the extrapolation. 



volumes agree within errors. Fig. [7] shows the agreement for a between the two volumes. As discussed in Section 
IIIIB1 there is an additional systematic error in calculating the nucleon mass on the ensemble with the lightest valence 
quark mass (am u = 0.005). For a conservative analysis we performed an extrapolation for a and P both with and 
without this lightest point. This gave a result which differed by 18% for a and 17% for (3 as shown in Fig. [5] We use 
this as an estimate of the error in extrapolating to the chiral limit. 

It should be noted that in our simulation, the strange quark mass is held fixed and hence in the extrapolation, 
only the light quarks are taken to the chiral limit. However, if we compare our result with the Nf = 2 result from 
|15| we see there is very good agreement (see Fig. For Nf = 2, the strange quark mass is effectively infinite, the 
agreement signifies that a and (3 have little dependence on the strange sea quark mass. 

For the NPR, we estimate a systematic error of 8% which is dominated by the error from truncating the pcrturbative 
expansion for the matching factor at order a 2 s in Eq. 1461 

Adding all of these uncertainties in quadrature, and together with the values for the matrix elements in Table [TTTT 
we estimate the low-energy parameters renormalised at [i = 2 GeV to be: 



a = -0.0112 ± 0.0012 (stat) ± 0.0022 (syst) GeV 3 (50) 
p = 0.0120 ± 0.0013( stat ) ± 0.0023( syst ) GeV 3 . (51) 

The results for various determinations of a are summarized in Fig. [5] The agreement between recent lattice 
computations suggests that lattice QCD is being successful at determining the low-energy constants describing nucleon 
decay with increasingly smaller systematic uncertainty. 

The indirect computation of the proton lifetime has a further non-linear systematic error, due to the use of chiral 
perturbation theory in a kinematic regime where the pion has a large momentum. The relevant matrix element has 
been computed using both the indirect and direct methods in Ref. [15| where sizeable differences were seen between 
the two methods. For the case of the matrix elements in Eq. [17] and 1181 the indirect method was found to give an 
estimate for the matrix element of about two times larger than the direct method. 

Finally, let us discuss one way to use our result to discriminate between GUTs. The proton partial decay width in 
a generic channel can be split into 



r = LEC 2 x Aq C d x ^1 G ut, 



(52) 
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FIG. 7: The LEC a measured on the two different volumes. There are no noticeable finite size effects. The chiral extrapolations 
from the two volumes are shown as white filled circles and also agree within errors. 
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FIG. 8: An extrapolation for a and /3 both with and without the value from the lightest valence quark mass point. This gives 
results differing by 18% for a and 17% for /3 



where LEC is the low-energy constant, a or f3, that we calculate earlier in this paper. Aqcd contains information 
from QCD parameters and Aqut contains all the information about the underlying high-energy theory, including 
constants from the GUT. For the Minimal SUSY SU(5) GUT, expressions for the lifetime have been calculated for 
several decay modes in Ref. [34| . 

Dimensional Analysis gives the value of the proton lifetime as ~ a GUT TO p/-^GUT- Taking Mqtjt ~ 10 16 [22j 
and obtaining aGUT by running the strong coupling up to the GUT scale gives Tn ~ 10~ 68 GeV. The natural scale 
for Agut is Mq^ t . The values of the low-energy constant that we have computed, together with the values of the 
quantities in Aqcd [22| and the experimental bounds on the proton life time (35l . [36| . allow us to put bounds on 
^4-gut for different decay modes as summarized in Table [VTl The bounds quoted are at a 68% confidence level. The 
different decay modes provide more or less stringent bounds on j4gut- The viability of any GUT can be assessed by 
checking wether the relevant bounds are satisfied. 
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FIG. 9: Summary of computations of the hadronic matrix element a, as given in Tab. [Vj Square points correspond to QCD 
model calculations, blue circles correspond to Nf — lattice QCD calculations, the green circle is from Nf = 2 and the result 
from our Nf = 2 + 1 calculation is shown in red. 



As a simple example, for the decay mode p 
Agut(p — > e + 7r°) is given by [13] to be 



via X boson exchange in the minimal SU(5) SUSY GUT, 



A 



GUT 



(p e+7r°) 



9jA 2 R 
M x 



\V ua 



(53) 



where 175 is the unified coupling at the GUT scale, Mx is the mass of the X boson w Maur, An is the renormalization 
factor and V u d a CKM matrix element. Using the value of Ar given in [37]|, we can put a bound on the X boson mass 
ofM x >5x 10 15 GeV. 

The decay widths of the channels involving colour triplet Higgs exchange can be calculated and involve the LEC /3 
(see Ref. The analysis in Ref. [l[ uses a conservative choice of j3 = 0.003GcV 3 at a scale of 1 GeV to constrain 

the mass of the colour triplet Higgs sufficiently to rule out the minimal SUSY SU(5) GUT. The higher value calculated 
in this work (running our value of (3 to a scale of 1 GeV gives (3 = 0.0109 ± 23 GeV 3 if we use Eq. Wf\ gives an even 
stronger constaint on the mass of the colour triplet Higgs and so confirms the fact that the minimal SUSY SU(5) 
GUT has been ruled out. 

The uncertainty on a 2 is 45% and on 1 is 43%. These are higher than the uncertainties on the factors Aqcd which 
for all channels is m 8%. A factor of sw 2 reduction in the uncertainty of the LECs would make them comparable with 
the uncertainties of Aqcd, which is dominated by the uncertainties of Z), F and f w . As Mx ~ \/a, an error of 45% 
on a 2 corresponds to an error of 11% on the bound for Mx- Reducing the uncertainty on a by a factor of two would 
therefore reduce the uncertainty from a on the bound for Mx to 6%. 
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|a| [GeV 3 ] 


\(3\ [GeV 3 ] 






Donoghue and Goldwich [4] 


0.003 




Bag model 




Thomas and McKellar [7] 


0.02 




Bag model 




Meljanac et al. [5] 


0.004 




Bag model 


QCD model 


loffe [2] 


0.009 




Sum rule 


calculation 


Krasnikov et al. [6] 


0.003 




Sum rule 




loffe and Smilga [8] 


0.006 




Sum rule 




Tomozawa [3] 


0.006 




Quark model 




Brodsky et al. [9] 


0.03 








Hara et al. [10] 


0.03 




WF, a = 0.11 fm 




Bowler et al. [11] 


0.013 


0.010 


WF, a = 0.22 fm 


Lattice QCD 


Gavela et al. [12] 


0.0056(8) 




WF, a = 0.09 fm 


N f =0 


JLQCD [13] 


0.015(1) 


0.014(1) 


WF, a = 0.09 fm 




CP-PACS & JLQCD [14] 


0.0090(09) (t^g) 


0.0096(09)(t^ ) 


WF, continuum limit 




Aoki et al. [15] 


0.0100(19) 


0.0108(21) 


DWF, a = 0.15 fm 


Lattice QCD 
N f =2 


Aoki et al. [15] 


0.0118(21) 


0.0118(21) 


DWF, a = 0.12 fm 


Lattice QCD 
Nf = 2 + 1 


This work 


0.0112(25) 


0.0120(26) 


DWF, a = 0.12 fm 



TABLE V: Comparison of the low energy parameter of the nucleon decay chiral Lagrangian a and (3 among various QCD 
model calculation, lattice results in the literatures and the results from this work. In lattice QCD calculations, WF and DWF 
mean Wilson and domain-wall fermions. The results for Nf = 2, and our results for Nf = 2 + 1 are shown with the total 
error consisting of statistical and systematic errors on the bare matrix element and renormalization constant. The errors on 
the results from Nf = are only statistical. 



Decay Mode 


Lifetime bound(yrs) 


Agvt bound (M GUT ) 


p — > e + 7r° 


> 8.2 x 10 33 


< 44 


p — > e + 7r° 


> 8.2 x 10 33 


< 37 


p -» K+u 


> 2.3 x 10 33 


< 76 


n -> K°D 


> 1.3 x 10 32 


< 733 



TABLE VI: The experimental proton partial lifetime bounds at 90% CL from [35|, H|| and the bound on Agut at a 68% CL that 
this lifetime bound implies. This bound is given in units of Mq^ t , the numbers quoted in the table are therefore dimensionless. 
The first line is for a decay mediated by a heavy gauge boson, the second and subsequent lines are for decays mediated by a 
colour triplet Higgs. 
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